clear
N=100;%number of k points

u1=0.0797*01;
u2=0.0975*01;
v=2.1354;
theta=1.09*pi/180; %acos(5953/5954);
% theta=acos(6337/6338);
D=0;
n=4;
q=8*pi/sqrt(3)*sin(theta/2);
Ga=4*pi/3*cos(theta/2)+q/2;
Ga=Ga*[1;0];

kx=q*linspace(-1,1,N)*0.8;
ky=kx;

B0=[0;0];Bx=[1;0];By=[0;1]; % magnetic field, unit is Tesla

%%
%from K to Gamma
for i=1:N
    tic
    for j=1:N
        i,j
        k=[kx(i);ky(j)];
    
        hp0=CM(u1,u2,v,theta,D,k-Ga,+1,n,B0);
        dp0=sort(eig(hp0));
        Ep0(j,i)=dp0(8*n^2);
        
        hn0=CM(u1,u2,v,theta,D,k+Ga,-1,n,B0);
        dn0=sort(eig(hn0));
        En0(j,i)=dn0(8*n^2);
        
        hpx=CM(u1,u2,v,theta,D,k-Ga,+1,n,Bx);
        dpx=sort(eig(hpx));
        Epx(j,i)=dpx(8*n^2);
        
        hnx=CM(u1,u2,v,theta,D,k+Ga,-1,n,Bx);
        dnx=sort(eig(hnx));
        Enx(j,i)=dnx(8*n^2);
        
        hpy=CM(u1,u2,v,theta,D,k-Ga,+1,n,By);
        dpy=sort(eig(hpy));
        Epy(j,i)=dpy(8*n^2);
        
        hny=CM(u1,u2,v,theta,D,k+Ga,-1,n,By);
        dny=sort(eig(hny));
        Eny(j,i)=dny(8*n^2);
        
    end
    toc
end
%%
save('FSHF.mat', 'kx', 'ky', 'Ep0', 'En0', 'Epx', 'Enx', 'Epy', 'Eny', 'q');

%%

load('FSHF.mat', 'kx', 'ky', 'Ep0', 'En0', 'Epx', 'Enx', 'Epy', 'Eny', 'q');

Kth = q/sqrt(3)/pi;

fsarea = @(c) polyarea(c(1,end-max(c(2,c(2,:)>=1))+1:end), c(2,end-max(c(2,c(2,:)>=1))+1:end));
hexarea = (Kth)^2 * 3*sqrt(3)/2;

theta1= 0:pi/3:2*pi;
figure(6997);
clf;
subplot(1,2,1)
hold on

Ef0 = (1.665-0.5352)/1000;
cp0=contour(kx/pi,ky/pi,Ep0,[Ef0, Ef0],'r','LineWidth',1);
fsarea(cp0)/hexarea

Efx = (1.665-0.5442)/1000;
cpx=contour(kx/pi,ky/pi,Epx,[Efx, Efx],'b','LineWidth',1);
fsarea(cpx)/hexarea

% contour(kx/pi,ky/pi,Epy,vv,'g','LineWidth',1)
colormap(white)
plot(q/sqrt(3)*sin(theta1)/pi,q/sqrt(3)*cos(theta1)/pi,'k--')
xlabel('k_x/\pi')
ylabel('k_y/\pi')
axis equal
title '+ valley'
subplot(1,2,2)
hold on

cn0=contour(kx/pi,ky/pi,En0,[Ef0, Ef0],'r','LineWidth',1);
fsarea(cn0)/hexarea

cnx=contour(kx/pi,ky/pi,Enx,[Efx, Efx],'b','LineWidth',1);
fsarea(cnx)/hexarea
% contour(kx/pi,ky/pi,Eny,vv,'g','LineWidth',1)
colormap(white)
% hold on
plot(q/sqrt(3)*sin(theta1)/pi,q/sqrt(3)*cos(theta1)/pi,'k--')
xlabel('k_x/\pi')
ylabel('k_y/\pi')
axis equal
title '- valley'

%% Analysis energy gap along FS

gv1 = zeros(N);
gv2 = gv1;
sigma = 0.0001;

for i=1:N
    for j=1:N
        e0 = En0(i,j);
        e5 = Ep0(i,j);
        e1 = Enx(i,j);
        e2 = Epx(end-i+1, end-j+1); % E(-k) of the opposite valley
        e3 = Enx(end-i+1, end-j+1);
        e4 = Epx(i,j); % E(-k) of the opposite valley
        
        gv1(i,j) = (e2-e1)*exp(-((e0-Ef0)/sigma)^2);
        gv2(i,j) = (e4-e3)*exp(-((e5-Ef0)/sigma)^2);
    end
end

gv = gv1 - gv2;

figure(6098);
clf;
hold on;
surf(kx/pi, ky/pi, zeros(N), gv1-gv2, 'EdgeColor', 'none'); axis tight;
plot(q/sqrt(3)*sin(theta1)/pi,q/sqrt(3)*cos(theta1)/pi,'k--')
% contour(kx/pi,ky/pi,En0,[Ef0, Ef0],'r','LineWidth',1);
colormap(blueredcmap());
axis equal

%% Export Data

fs = @(c) c(:,end-max(c(2,c(2,:)>=1))+1:end)/Kth;

fp = fopen('FS1p09HF.txt', 'w');
fprintf(fp, '# K1\n');
fprintf(fp, '%g\t%g\n', fs(cp0));
fprintf(fp, '\n\n');

fprintf(fp, '# K2\n');
fprintf(fp, '%g\t%g\n', fs(cn0));
fprintf(fp, '\n\n');

fprintf(fp, '# K1x\n');
fprintf(fp, '%g\t%g\n', fs(cpx));
fprintf(fp, '\n\n');

fprintf(fp, '# K2x\n');
fprintf(fp, '%g\t%g\n', fs(cnx));
fprintf(fp, '\n\n');

fprintf(fp, '# GV\n');
for i=1:N
    fprintf(fp, '%g\t%g\t%g\n', [repmat(kx(i),N,1)/Kth/pi, ky(:)/Kth/pi, gv1(i,:)']');
    fprintf(fp, '\n');
end
fprintf(fp, '\n');

fclose(fp);